home *** CD-ROM | disk | FTP | other *** search
/ Fritz: All Fritz / All Fritz.zip / All Fritz / FILES / PROGMISC / PCSSP.LZH / PC-SSP.ZIP / POLYROOT.ZIP / PQFB.FOR < prev    next >
Text File  |  1985-11-29  |  7KB  |  238 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE PQFB
  5. C
  6. C        PURPOSE
  7. C           TO FIND AN APPROXIMATION Q(X)=Q1+Q2*X+X*X TO A QUADRATIC
  8. C           FACTOR OF A GIVEN POLYNOMIAL P(X) WITH REAL COEFFICIENTS.
  9. C
  10. C        USAGE
  11. C           CALL PQFB(C,IC,Q,LIM,IER)
  12. C
  13. C        DESCRIPTION OF PARAMETERS
  14. C           C   - INPUT VECTOR CONTAINING THE COEFFICIENTS OF P(X) -
  15. C                 C(1) IS THE CONSTANT TERM (DIMENSION IC)
  16. C           IC  - DIMENSION OF C
  17. C           Q   - VECTOR OF DIMENSION 4 - ON INPUT Q(1) AND Q(2) MUST
  18. C                 CONTAIN INITIAL GUESSES FOR Q1 AND Q2 - ON RETURN Q(1)
  19. C                 AND Q(2) CONTAIN THE REFINED COEFFICIENTS Q1 AND Q2 OF
  20. C                 Q(X), WHILE Q(3) AND Q(4) CONTAIN THE COEFFICIENTS A
  21. C                 AND B OF A+B*X, WHICH IS THE REMAINDER OF THE QUOTIENT
  22. C                 OF P(X) BY Q(X)
  23. C           LIM - INPUT VALUE SPECIFYING THE MAXIMUM NUMBER OF
  24. C                 ITERATIONS TO BE PERFORMED
  25. C           IER - RESULTING ERROR PARAMETER (SEE REMARKS)
  26. C                 IER= 0 - NO ERROR
  27. C                 IER= 1 - NO CONVERGENCE WITHIN LIM ITERATIONS
  28. C                 IER=-1 - THE POLYNOMIAL P(X) IS CONSTANT OR UNDEFINED
  29. C                          - OR OVERFLOW OCCURRED IN NORMALIZING P(X)
  30. C                 IER=-2 - THE POLYNOMIAL P(X) IS OF DEGREE 1
  31. C                 IER=-3 - NO FURTHER REFINEMENT OF THE APPROXIMATION TO
  32. C                          A QUADRATIC FACTOR IS FEASIBLE, DUE TO EITHER
  33. C                          DIVISION BY 0, OVERFLOW OR AN INITIAL GUESS
  34. C                          THAT IS NOT SUFFICIENTLY CLOSE TO A FACTOR OF
  35. C                          P(X)
  36. C
  37. C        REMARKS
  38. C           (1)  IF IER=-1 THERE IS NO COMPUTATION OTHER THAN THE
  39. C                POSSIBLE NORMALIZATION OF C.
  40. C           (2)  IF IER=-2 THERE IS NO COMPUTATION OTHER THAN THE
  41. C                NORMALIZATION OF C.
  42. C           (3)  IF IER =-3  IT IS SUGGESTED THAT A NEW INITIAL GUESS BE
  43. C                MADE FOR A QUADRATIC FACTOR.  Q, HOWEVER, WILL CONTAIN
  44. C                THE VALUES ASSOCIATED WITH THE ITERATION THAT YIELDED
  45. C                THE SMALLEST NORM OF THE MODIFIED LINEAR REMAINDER.
  46. C           (4)  IF IER=1, THEN, ALTHOUGH THE NUMBER OF ITERATIONS LIM
  47. C                WAS TOO SMALL TO INDICATE CONVERGENCE, NO OTHER PROB-
  48. C                LEMS HAVE BEEN DETECTED, AND Q WILL CONTAIN THE VALUES
  49. C                ASSOCIATED WITH THE ITERATION THAT YIELDED THE SMALLEST
  50. C                NORM OF THE MODIFIED LINEAR REMAINDER.
  51. C           (5)  FOR COMPLETE DETAIL SEE THE DOCUMENTATION FOR
  52. C                SUBROUTINES PQFB AND DPQFB.
  53. C
  54. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  55. C           NONE
  56. C
  57. C        METHOD
  58. C           COMPUTATION IS BASED ON BAIRSTOW'S ITERATIVE METHOD.  (SEE
  59. C           WILKINSON, J.H., THE EVALUATION OF THE ZEROS OF ILL-CON-
  60. C           DITIONED POLYNOMIALS (PART ONE AND TWO), NUMERISCHE MATHE-
  61. C           MATIK, VOL.1 (1959), PP. 150-180, OR HILDEBRAND, F.B.,
  62. C           INTRODUCTION TO NUMERICAL ANALYSIS, MC GRAW-HILL, NEW YORK/
  63. C           TORONTO/LONDON, 1956, PP. 472-476.)
  64. C
  65. C     ..................................................................
  66. C
  67.       SUBROUTINE PQFB(C,IC,Q,LIM,IER)
  68. C
  69. C
  70.       DIMENSION C(1),Q(1)
  71. C
  72. C        TEST ON LEADING ZERO COEFFICIENTS
  73.       IER=0
  74.       J=IC+1
  75.     1 J=J-1
  76.       IF(J-1)40,40,2
  77.     2 IF(C(J))3,1,3
  78. C
  79. C        NORMALIZATION OF REMAINING COEFFICIENTS
  80.     3 A=C(J)
  81.       IF(A-1.)4,6,4
  82.     4 DO 5 I=1,J
  83.       C(I)=C(I)/A
  84.       CALL OVERFL(N)
  85.       IF(N-2)40,5,5
  86.     5 CONTINUE
  87. C
  88. C        TEST ON NECESSITY OF BAIRSTOW ITERATION
  89.     6 IF(J-3)41,38,7
  90. C
  91. C        PREPARE BAIRSTOW ITERATION
  92.     7 EPS=1.E-6
  93.       EPS1=1.E-3
  94.       L=0
  95.       LL=0
  96.       Q1=Q(1)
  97.       Q2=Q(2)
  98.       QQ1=0.
  99.       QQ2=0.
  100.       AA=C(1)
  101.       BB=C(2)
  102.       CB=ABS(AA)
  103.       CA=ABS(BB)
  104.       IF(CB-CA)8,9,10
  105.     8 CC=CB+CB
  106.       CB=CB/CA
  107.       CA=1.
  108.       GO TO 11
  109.     9 CC=CA+CA
  110.       CA=1.
  111.       CB=1.
  112.       GO TO 11
  113.    10 CC=CA+CA
  114.       CA=CA/CB
  115.       CB=1.
  116.    11 CD=CC*.1
  117. C
  118. C        START BAIRSTOW ITERATION
  119. C        PREPARE NESTED MULTIPLICATION
  120.    12 A=0.
  121.       B=A
  122.       A1=A
  123.       B1=A
  124.       I=J
  125.       QQQ1=Q1
  126.       QQQ2=Q2
  127.       DQ1=HH
  128.       DQ2=H
  129. C
  130. C        START NESTED MULTIPLICATION
  131.    13 H=-Q1*B-Q2*A+C(I)
  132.       CALL OVERFL(N)
  133.       IF(N-2)42,14,14
  134.    14 B=A
  135.       A=H
  136.       I=I-1
  137.       IF(I-1)18,15,16
  138.    15 H=0.
  139.    16 H=-Q1*B1-Q2*A1+H
  140.       CALL OVERFL(N)
  141.       IF(N-2)42,17,17
  142.    17 C1=B1
  143.       B1=A1
  144.       A1=H
  145.       GO TO 13
  146. C        END OF NESTED MULTIPLICATION
  147. C
  148. C        TEST ON SATISFACTORY ACCURACY
  149.    18 H=CA*ABS(A)+CB*ABS(B)
  150.       IF(LL)19,19,39
  151.    19 L=L+1
  152.       IF(ABS(A)-EPS*ABS(C(1)))20,20,21
  153.    20 IF(ABS(B)-EPS*ABS(C(2)))39,39,21
  154. C
  155. C        TEST ON LINEAR REMAINDER OF MINIMUM NORM
  156.    21 IF(H-CC)22,22,23
  157.    22 AA=A
  158.       BB=B
  159.       CC=H
  160.       QQ1=Q1
  161.       QQ2=Q2
  162. C
  163. C        TEST ON LAST ITERATION STEP
  164.    23 IF(L-LIM)28,28,24
  165. C
  166. C        TEST ON RESTART OF BAIRSTOW ITERATION WITH ZERO INITIAL GUESS
  167.    24 IF(H-CD)43,43,25
  168.    25 IF(Q(1))27,26,27
  169.    26 IF(Q(2))27,42,27
  170.    27 Q(1)=0.
  171.       Q(2)=0.
  172.       GO TO 7
  173. C
  174. C        PERFORM ITERATION STEP
  175.    28 HH=AMAX1(ABS(A1),ABS(B1),ABS(C1))
  176.       IF(HH)42,42,29
  177.    29 A1=A1/HH
  178.       B1=B1/HH
  179.       C1=C1/HH
  180.       H=A1*C1-B1*B1
  181.       IF(H)30,42,30
  182.    30 A=A/HH
  183.       B=B/HH
  184.       HH=(B*A1-A*B1)/H
  185.       H=(A*C1-B*B1)/H
  186.       Q1=Q1+HH
  187.       Q2=Q2+H
  188. C        END OF ITERATION STEP
  189. C
  190. C        TEST ON SATISFACTORY RELATIVE ERROR OF ITERATED VALUES
  191.       IF(ABS(HH)-EPS*ABS(Q1))31,31,33
  192.    31 IF(ABS(H)-EPS*ABS(Q2))32,32,33
  193.    32 LL=1
  194.       GO TO 12
  195. C
  196. C        TEST ON DECREASING RELATIVE ERRORS
  197.    33 IF(L-1)12,12,34
  198.    34 IF(ABS(HH)-EPS1*ABS(Q1))35,35,12
  199.    35 IF(ABS(H)-EPS1*ABS(Q2))36,36,12
  200.    36 IF(ABS(QQQ1*HH)-ABS(Q1*DQ1))37,44,44
  201.    37 IF(ABS(QQQ2*H)-ABS(Q2*DQ2))12,44,44
  202. C        END OF BAIRSTOW ITERATION
  203. C
  204. C        EXIT IN CASE OF QUADRATIC POLYNOMIAL
  205.    38 Q(1)=C(1)
  206.       Q(2)=C(2)
  207.       Q(3)=0.
  208.       Q(4)=0.
  209.       RETURN
  210. C
  211. C        EXIT IN CASE OF SUFFICIENT ACCURACY
  212.    39 Q(1)=Q1
  213.       Q(2)=Q2
  214.       Q(3)=A
  215.       Q(4)=B
  216.       RETURN
  217. C
  218. C        ERROR EXIT IN CASE OF ZERO OR CONSTANT POLYNOMIAL
  219.    40 IER=-1
  220.       RETURN
  221. C
  222. C        ERROR EXIT IN CASE OF LINEAR POLYNOMIAL
  223.    41 IER=-2
  224.       RETURN
  225. C
  226. C        ERROR EXIT IN CASE OF NONREFINED QUADRATIC FACTOR
  227.    42 IER=-3
  228.       GO TO 44
  229. C
  230. C        ERROR EXIT IN CASE OF UNSATISFACTORY ACCURACY
  231.    43 IER=1
  232.    44 Q(1)=QQ1
  233.       Q(2)=QQ2
  234.       Q(3)=AA
  235.       Q(4)=BB
  236.       RETURN
  237.       END
  238.